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ABSTRACT 

We report on the discovery of a planetary system with a close-in transiting hot Jupiter on a near 
circular orbit and a massive outer planet on a highly eccentric orbit. The inner planet, HAT-P-13b, 
transits the bright V=10.622 G4 dwarf star GSC 3416-00543 every P = 2.916260±0.000010 days, with 
transit epoch T c = 2454779.92979 ± 0.00038 (BJD) and duration 0.1345 ± 0.0017 d. The outer planet, 
HAT-P-13c orbits the star with P2 — 428.5±3.0 days and nominal transit center (assuming zero impact 
parameter) of T 2c = 2454870.4 ± 1.8 (BJD) or time of periastron passage T 2iPe „ = 2454890.05 ± 0.48 
(BJD). Transits of the outer planet have not been observed, and may not be present. The host star 
has a mass of L22~t ) 'j ) u ' Mq, radius of 1.56 ± 0.08 Rq, effective temperature 5653 ± 90 K, and is rather 
metal rich with [Fe/H] = +0.41 ± 0.08. The inner planetary companion has a mass of 0.853^0.046 -^J; 
and radius of 1.281 ±0.079 Rj yielding a mean density of 0.498loo69 § cm -3 . The outer companion has 
ui2 sini2 = 15.2 ± 1.0 Mj, and orbits on a highly eccentric orbit of e 2 = 0.691 ± 0.018. While we have 
not detected significant transit timing variations of HAT-P-13b, due to gravitational and light-travel 
time effects, future observations will constrain the orbital inclination of HAT-P-13c, along with its 
mutual inclination to HAT-P-13b. The HAT-P-13 (b,c) double-planet system may prove extremely 
valuable for theoretical studies of the formation and dynamics of planetary systems. 
Subject headings: planetary systems — stars: individual (HAT-P-13, GSC 3416-00543) techniques: 
spectroscopic, photometric 
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1. INTRODUCTION 

Radial velocity (RV) surveys have shown that multiple- 
planet stellar systems are common. For example, Wright 
et al. (2007) concluded that the occurrence of additional 
planets among stars already having one known planet 
must be greater than 30%. Thus, one would expect that 
out of the ~ 50 published transiting extrasolar planet 
(TEP) systems 10 , there should be a number of systems 
with additional companions; these companions should 
make their presence known through radial velocity vari- 
ations of the parent stars, even if they do not them- 
selves transit. The fact that so far no TEPs in multiple 
planet systems have been reported is somewhat surpris- 
ing based on the statistics from RV surveys (see also Fab- 
rycky 2009). Recently Smith et al. (2009) searched the 
light curves of 24 transiting planets for outer compan- 
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ions, finding no evidence for a double transiting system. 

The Hungarian-made Automated Telescope Network 
(HATNet) survey (Bakos et al. 2002, 2004) has been 
a major contributor to the discovery of TEPs. Oper- 
ational since 2003, it has covered approximately 10% of 
the Northern sky, searching for TEPs around bright stars 
(8 ^ I ;$ 12.5 mag). HATNet operates six wide field in- 
struments: four at the Fred Lawrence Whipple Obser- 
vatory (FLWO) in Arizona, and two on the roof of the 
Submillimeter Array hangar (SMA) of SAO at Hawaii. 
Since 2006, HATNet has announced and published 12 
TEPs. A study similar to that of Smith et al. (2009) has 
been carried out on 9 known transiting planets from the 
HATNet project by Fabrycky (private communication) 
with a similar null result. 

In this work we report on the 13th discovery of HAT- 
Net, the detection of the first known system with a tran- 
siting inner planet (HAT-P-13b) which also contains a 
second, outer planet (HAT-P-13c), as detected by the 
radial velocity variation of the host star. There have 
been examples of transiting systems where the RV vari- 
ations do show a long term trend, such as HAT-P-7b 
(Winn et al. 2009) and HAT-P-llb (Bakos et al. 2009), 
but no orbital solution has yet been presented for any 
of these outer companions, simply because there has 
not been a long enough timespan to cover the period, 
or at least observe the long term trend changing sign. 
While no transits of HAT-P-13c have yet been detected, 
the probability that the companion actually transits the 
star is non-negligible if the orbits of the two planets are 
nearly co-planar (the pure geometric transit probability 
for HAT-P-13c is 1.3%, see § 4.3). The system is par- 
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ticularly interesting because the outer planet has both a 
high eccentricity and a very high mass. These properties, 
in turn, should induce transit timing variations (TTVs) 
of the inner planet of the order of 10 seconds (standard 
deviation, Agol et al. 2005). Such TTVs may be used 
to constrain the orbital parameters of the outer planet, 
including the inclination with respect to the line of sight 
and with respect to the orbital plane of the inner planet. 

In § 2 of this paper we summarize the observations, 
including the photometric detection, and follow-up ob- 
servations. § 3 describes the analysis of the data, such as 
the stellar parameter determination (§ 3.1), blend mod- 
eling (§ 3.2), and global modeling of the data (§ 3.3). We 
discuss our findings in § 4. 

2. OBSERVATIONS 

2.1. Photometric detection 

The transits of HAT-P-13b were detected with the 
HAT- 5 telescope. The region around GSC 3416-00543, a 
field internally labeled as 136, was observed on a nightly 
basis between 2005 November 25 and 2006 May 20, 
whenever weather conditions permitted. We gathered 
4021 exposures of 5 minutes at a 5.5-minute cadence. 
Each image contained approximately 20000 stars down to 
I ~ 14.0. For the brightest stars in the field we achieved 
a per-image photometric precision of 3.1 mmag. 
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Fig. 1. — The unbinned light curve of HAT-P-13 including all 
4021 instrumental / band measurements obtained with the HAT-5 
telescope of HATNet (see text for details), and folded with the pe- 
riod of P = 2.9162595 days (which is the result of the fit described 
in § 3). 

The calibration of the HATNet frames was done uti- 
lizing standard procedures. The calibrated frames were 
then subjected to star detection and astrometry, as de- 
scribed in Pal & Bakos (2006). Aperture photometry 
was performed on each image at the stellar centroids de- 
rived from the 2MASS catalog (Skrutskie et al. 2006) 
and the individual astrometrical solutions. Description 
of numerous details related to the reduction are also 
given in Pal (2009). The resulting light curves were 
decorrelated against trends using the External Parame- 
ter Decorrelation technique in "constant" mode (EPD, 
see Bakos et al. 2009) and the Trend Filtering Algorithm 
(TFA, see Kovacs et al. 2005). The light curves were 
searched for periodic box-like signals using the Box Least 
Squares method (BLS, see Kovacs et al. 2002). We de- 
tected a significant signal in the light curve of GSC 3416- 
00543 (also known as 2MASS 08393180+4721073; a = 
08 h 39 m 31.82s, S = +47°21'07.4"; J2000), with a depth 
of ~ 6.2 mmag, and a period of P — 2.9163 days. The dip 
significance parameter (Kovacs et al. 2002) of the signal 
was DSP = 17. The dip had a relative duration (first to 
last contact) of q « 0.0461 ± 0.0006, corresponding to a 
total duration of Pq « 3.228 ± 0.040 hours (see Fig. 1). 



2.2. Reconnaissance Spectroscopy 

One of the important tools used in our survey for estab- 
lishing whether the transit-feature in the light curve of a 
candidate is due to astrophysical phenomena other than 
a planet transiting a star is the CfA Digital Speedome- 
ter (DS; Latham 1992), mounted on the FLWO 1.5m 
telescope. This yields high-resolution spectra with low 
signal-to-noise ratio sufficient to derive radial velocities 
with moderate precision (roughly lkms -1 ), and to de- 
termine the effective temperature and surface gravity of 
the host star. With this facility we are able to weed 
out certain false alarms, such as eclipsing binaries and 
multiple stellar systems. 

We obtained 7 spectra spanning 37 days with the DS. 
The RV measurements of HAT-P-13 showed an rms resid- 
ual of 0.68 km s _1 , consistent with no detectable RV vari- 
ation. The spectra were single- lined, showing no evidence 
for more than one star in the system. Atmospheric pa- 
rameters for the star, including the initial estimates of 
effective temperature T e s* = 5500 ± 100 K, surface grav- 
ity log = 4.0 ±0.25 (log cgs), and projected rotational 
velocity wsini = 0.5 ± l.Okms -1 , were derived as de- 
scribed by Torres et al. (2002). The effective tempera- 
ture and surface gravity correspond to a G4 dwarf (Cox 
2000). The mean line-of-sight velocity of HAT-P-13 is 
+ 14.69 ± 0.68 kins" 1 . 

2.3. High resolution, high S/N spectroscopy and the 
search for radial velocity signal components 

Given the significant transit detection by HATNet, 
and the encouraging DS results, we proceeded with the 
follow-up of this candidate by obtaining high-resolution 
and high S/N spectra to characterize the radial velocity 
variations and to determine the stellar parameters with 
higher precision. Using the HIRES instrument (Vogt et 
al. 1994) on the Keck I telescope located on Mauna Kea, 
Hawaii, we obtained 30 exposures with an iodine cell, 
plus three iodine-free templates. The first template had 
low signal-to-noise ratio, thus we repeated the template 
observations during a later run, and acquired two high 
quality templates. We used the last two templates in 
the analysis. The observations were made between 2008 
March 22 and 2009 June 5. The relative radial velocity 
(RV) measurements are listed in Tab. 1, and shown in 
Fig. 2. 

Observations and reductions have been carried out in 
an identical way to that described in earlier HATNet 
discovery papers, such as Bakos et al. (2009). References 
for the Keck iodine cell observations and the reduction of 
the radial velocities are given in Marcy & Butler (1992) 
and Butler et al. (1996). 

Initial fits of these RVs to a single-planet Keplcrian 
orbit were quite satisfactory but soon revealed a slight 
residual trend that became more significant with time. 
This is the reason for the observations extending over 
more than one year, as opposed to just a few months as 
necessary to confirm simpler purely sinusoidal variations 
seen in other transiting planets discovered by HATNet. 
Eventually we noticed that the residual trend reversed 
(Fig. 2, top), a clear sign of a coherent motion most likely 
due to a more distant body in the same system, possi- 
bly a massive planet. Preliminary two-planet orbital so- 
lutions provided a much improved fit (although with a 
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TABLE 1 

Relative radial velocity, bisector and activity index measurements of 

HAT-P-13. 
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0.0042 
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929.84447 
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2.88 


-19.31 
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0.0040 
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955.86964 


-90.43 
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-13.92 
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0.0040 


0.00003 
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95.55 
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-5.32 
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0.0039 


0.00003 


963.85163 


-20.24 


1.62 


-5.27 


7.21 


0.0041 
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1.47 


4.02 
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0.0040 
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0.0041 
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5.66 


0.0040 
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6.37 


5.88 


0.0042 


0.00002 


986.76358 
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1.73 


3.95 


6.31 


0.0042 


0.00002 
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a The fitted zero-point that is on an arbitrary scale (denoted as 7 re ; in Tab. 4) 
has been subtracted from the velocities. 
b The values for cr RV do not include the jitter. 
c S values are on a relative scale. 



weakly determined outer period due to the short dura- 
tion of the observations), and more importantly, held up 
after additional observations. With the data so far gath- 
ered and presented in this work, a false alarm probabil- 
ity (FAP) analysis for the Keplerian orbit of the second 
planet yielded an FAP of 0.00001. 

A more thorough analysis using a two-component least 
squares Fourier fit (see Barning 1963, for the single- 
component case), with one component fixed to the known 
frequency of the short-period inner planet, also confirmed 
the existence of the long period component, and indi- 
cated that the latter is due to a highly non-sinusoidal 
motion. A number of other low frequency peaks were 
eliminated as being aliases yielding worse quality fits. A 
three-harmonic representation of the fit yielded an RMS 
of 5ms" 1 , fairly close to the value expected from the 
formal errors. Full modeling of this complex motion 
of two Keplerian orbits superposed, simultaneously with 
the photometry, is described in our global fit of § 3.3. 

2.4. Photometric follow-up observations 

We observed 7 transit events of HAT-P-13 with the 
KeplerCam CCD on the FLWO 1.2 m telescope between 
2008 April 24 and 2009 May 8. All observations were 
carried out in the i band, and the typical exposure time 
was 15 seconds. The reduction of the images, including 



calibration, astrometry, and photometry, was performed 
as described in Bakos et al. (2009). 

We performed EPD and TFA against trends simulta- 
neously with the light curve modeling (for more details, 
see § 3). From the series of apertures, for each night, 
we chose the one yielding the smallest residual rms for 
producing the final light curve. These are shown in the 
upper plot of Fig. 3, superimposed with our best fit tran- 
sit light curve models (see also § 3). 

3. ANALYSIS 
3.1. Properties of the parent star 

We derived the initial stellar atmospheric parameters 
by using the first template spectrum obtained with the 
Keck/HIRES instrument. We used the SME package of 
Valenti & Piskunov (1996) along with the atomic-line 
database of Valenti & Fischer (2005), which yielded the 
following initial values and uncertainties (which we have 
conservatively increased, to include our estimates of the 
systematic errors): effective temperature Toff* = 5760 ± 
90K, stellar surface gravity logg* = 4.28 ± 0.10 (cgs), 
mctallicity [Fe/H] = +0.50 ± 0.08 dex, and projected ro- 
tational velocity vsmi = 1.9 ± 0.5kms _1 . 

Further analysis of the spectra has shown that the host 
star is chromospherically quiet with \ogR' HK = —5.10 
and S — 0.14 (on an absolute scale). 
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Fig. 2. — (Top) Radial- velocity measurements from Keck for 
HAT-P-13, along with the best 2-planet orbital fit, shown as a 
function of BJD (see § 3). The center-of-mass velocity has been 
subtracted. (Second panel) Phased residuals after subtracting the 
orbital fit (also see § 3). The rms variation of the residuals is about 
3.4ms , requiring a jitter of 3.0ms" 1 to be added in quadrature 
to the individual errors to yield a r educed \ 2 of 1.0. The error-bars 
in this panel have been inflated accordingly. Note the different ver- 
tical scale of the panels. (Third panel) Orbit of the outer planet 
after subtracting the orbital fit of the inner planet from the data. 
Zero phase is defined by the fictitious transit midpoint of the sec- 
ond planet (denoted as T c 2, where c is the common subscript for 
"center", and "2" refers to the second planet). The error- bars (in- 
flated by the jitter) are smaller than the size of the points. (Fourth 
panel) Orbit of the inner planet after subtracting the orbital fit of 
the outer planet. Zero-phase is defined by the transit midpoint of 
HAT-P-13b (denoted as T c i). The error-bars are smaller than the 
size of the points. (Fifth panel) Bisector spans (BS) for the Keck 
spectra including the two template spectra (§ 3.2). The mean value 
has been subtracted. (Bottom) Relative S activity index values for 
the Keck spectra (see Hartman et al. 2009). Open circles in the 
phase-plots indicate data that appears twice due to plotting outside 
the [0,1] phase domain. 
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Fig. 3. — Unbinncd instrumental i band transit light curves, 
acquired with KeplerCam at the FLWO 1.2 m telescope on seven 
different dates. If the first full transit is assigned Ntr = (2008 
November 8/9 MST, the third event from the top), then the follow- 
up light curves have Ntr = -68, -1, 0, 1, 24, 35, 62 from top to 
bottom. Superimposed are the best-fit transit model light curves as 
described in § 3. In the bottom of the figure we show the residuals 
from the fit. Error-bars represent the photon and background shot- 
noise, plus the readout noise. 



TABLE 2 

Photometry follow-up of HAT-P-13 



BJD 


Mag a 


"Mag 


Mag(orig) Filter 


(2,400,000+) 






54581.66038 


0.00550 


0.00080 


0.03620 i 


54581.66070 


0.00763 


0.00080 


0.03652 i 


54581.66104 


0.00877 


0.00080 


0.03686 i 


54581.66138 


0.00725 


0.00080 


0.03720 i 


54581.66171 


0.00670 


0.00080 


0.03753 i 


54581.66205 


0.01004 


0.00080 


0.03787 i 


54581.66237 


0.00673 


0.00080 


0.03819 i 


54581.66271 


0.00787 


0.00080 


0.03853 i 



Note. — This table is presented in its entirety in the elec- 
tronic edition of the Astrophysical Journal. A portion is shown 
here for guidance regarding its form and content. 

a The out-of-transit level has been subtracted. These magni- 
tudes have been subjected to the EPD and TEA procedures (in 
"ELTG" mode), carried out simultaneously with the transit fit. 

To determine the stellar properties via a set of 
isochrones, we used three parameters: the stellar effective 
temperature, the metallicity, and the normalized semi- 
major axis a/ R* (or related mean stellar density p*). We 
note that another possible parameter in place of a/i?* 
would be the stellar surface gravity, but in the case of 
planetary transits the a/ R± quantity typically imposes a 
stronger constraint on the stellar models (Sozzetti et al. 
2007). (The validity of this assumption, namely that the 
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adequate physical model describing our data is a plane- 
tary transit, as opposed to e.g. a blend, is shown later in 
§ 3.2.) With quadratic limb darkening coefficients (listed 
in Tab. 3) from Claret (2004) appropriate for the values 
of T e ff*, [Fe/H], and loggias derived from the SME anal- 
ysis, we performed a global modeling of the data (§ 3.3), 
yielding a full Monte-Carlo distribution of a/i?+. For 
Toff* and [Fe/H]we adopted normal distributions in the 
Monte Carlo analysis, with dispersions equal to the l-a 
uncertainties from the initial SME analysis. 

For each combination within the large (~ 20000) set 
of a/Ri,, Teff*, [Fe/H] values, we searched the stellar 
isochrones of the Yi et al. (2001) models for the best fit 
stellar model parameters (M*, R*, age, logg*, etc.). We 
derived the mean values and uncertainties of the physical 
parameters based on their a posteriori distribution (see 
e.g. Pal et al. 2008). 

We then repeated the SME analysis by fixing the stellar 
surface gravity to the refined value of logg* = 4.13±0.04 
based on the isochrone search, and only adjusting T c ff*, 
[Fe/H] and vsini. The SME analysis was performed on 
the first, weaker template observation, and also on the 
second and third, higher S/N pair of template observa- 
tions taken by Keck/HIRES. While the pair of high S/N 
templates were acquired very close in time, and the re- 
spective SME values were consistent to within a small 
fraction of the formal error-bars, they were also consis- 
tent to within l-a with the values based on the weaker 
template that was taken much earlier. This consistency 
reassured that our stellar atmospheric parameter deter- 
mination is robust, and the error-bars are realistic. Be- 
cause the second two templates were of better quality, 
we adopted the SME values found from these spectra 
with simple averaging, yielding T g-* = 5653 ± 90 K, 
[Fe/H] = +0.41 ± 0.08 and vsini = 2.9 ± 1.0 km s" 1 . 
We adopted these as the final atmospheric parameters 
for the star. We then also repeated the isochrone search 
for stellar parameters, obtaining M+=l. 219^0 099 -^©> 
#*=1.559± 0.082 R Q and L*=2.22± 0.31 L Q . These are 
summarized in Tab. 3, along with other stellar proper- 
ties. Model isochrones from Yi et al. (2001) for metal- 
licity [Fe/H] =+0.41 are plotted in Fig. 4, with the final 
choice of effective temperature T e fj* and a/R* marked, 
and encircled by the l-a and 2-a confidence ellipsoids. 
The second SME iteration at fixed stellar surface gravity 
(as determined from a/R+) changed the metallicity and 
stellar temperature in such a way that the new (T e ff*, 
a/R±) values now fall in a more complex region of the 
isochrones, as compared to the initial SME values, allow- 
ing for multiple solutions (see Fig. 4, where the original 
SME values are marked with a triangle). 11 The distribu- 
tion of stellar age becomes bimodal with the dominant 
peak in the histogram (not shown) being at 5.0 Gyr, and 
a smaller peak (by a factor of 5) at 7.3 Gyr. This corre- 
sponds to a slightly bimodal mass distribution with the 
dominant peak at ~ 1.23 M©, and much smaller peak 
around ~ 1.13 M©. The asymmetric error-bars given 
in Tab. 3 for the mass and age account for the double- 
peaked distribution. 

The stellar evolution modeling also yields the abso- 

11 The reason for the intersecting isochrones is the "kink" on the 
evolutionary tracks for > 1.0 Mq stars evolving off the main 
sequence. 
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Fig. 4. — Stellar isochrones from Yi et al. (2001) for metallicity 
[Fe/H]=+0.41 and ages between 3.6 and 8.0 Gyr in steps of 0.4 Gyr. 
The final choice of T c ff* and a/R t are marked and encircled by 
the l-a and 2-cr confidence ellipsoids. The open triangle denotes 
the Teff+ja/R* point found during the first SME iteration, before 
refining the stellar surface gravity. 

lute magnitudes and colors in various photometric pass- 
bands. We used the apparent magnitudes from the 
2MASS catalogue (Skrutskic et al. 2006) to determine 
the distance of the system, after conversion to the ESO 
system of the models. The reported magnitudes for this 
star are J 2 mass = 9.328 + 0.018, H 2 mass = 9.058 + 0.017 
and i^2MASs = 8.975 ± 0.017, which we transformed to 
J ESO = 9.396 + 0.022, H ESO = 9.070 + 0.021 and K ESO = 
9.018+0.018, following Carpenter (see 2001). These yield 
a color of (J — K) = 0.378 ± 0.028, fully consistent with 
the expected, isochrone-based (J — K)i so = 0.374 + 0.015. 
We thus relied on the 2MASS K apparent magnitude and 
the Mk = 2.36 ± 0.12 absolute magnitude derived from 
the above-mentioned modeling to determine the distance: 
214 ± 12 pc, assuming no reddening. The K band was 
chosen because it is the longest wavelength band-pass 
with the smallest expected discrepancies due to molecu- 
lar lines in the spectrum of the star. 

3.2. Excluding blend scenarios 

Following Torres et al. (2007), we explored the pos- 
sibility that the measured radial velocities arc not due 
to the (multiple) planet-induced orbital motion of the 
star, but are instead caused by distortions in the spec- 
tral line profiles. This could be due to contamination 
from a nearby unresolved eclipsing binary, in this case 
presumably with a second companion producing the RV 
signal corresponding to HAT-P-13c. A bisector analysis 
based on the Keck spectra was done as described in §3 
of Hartman et al. (2009). 

We detect no significant variation in the bisector spans 
(see Fig. 2, fifth panel). The correlation between the ra- 
dial velocity and the bisector variations is also insignifi- 
cant. Therefore, we conclude that the velocity variations 
of the host star are real, and can be interpreted as be- 
ing due to a close-in planet, with the added complication 
from an outer object that we account for in the modeling 
that follows. Because of the negligible bisector variations 
that show no correlation with the radial velocities, we 
found no need to perform detailed blend modeling of hi- 
erarchical triple (quadruple) scenarios, such as that done 
for the case of HAT-P-12b (Hartman et al. 2009). 

3.3. Global modeling of the data 
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TABLE 3 
Stellar parameters for HAT-P-13 



Parameter Value Source 



T eff t (K) 5653 ± 90 SME a 

[Fe/H] (dex) . . . +0.41 ± 0.08 SME 

vsini (bus- 1 ) 2.9 ±1.0 SME 

fmac (kms -1 ) . 3.83 SME 

i) mic (kms- 1 ). . 0.85 SME 

7rv (kms- 1 ) . . +14.69±0.68 DS 

7i (i) 0.3060 SME+Claret 

7 2{i) 0.3229 SME+Claret 

M* (Mq) 1.22±° °g YY+a/R*+SME c 

R* (Rq) 1.56 ±0.08 YY+a/R*+SME 

logs* (cgs) .... 4.13 ±0.04 YY+a/R*+SME 

L* (L Q ) 2.22 ±0.31 YY+a/R*+SME 

V (mag) 10.622 TASS d 

B-V (mag) . . 0.73 ± 0.03 YY+«/-R*+SME 

My (mag) 3.97 ±0.17 YY+a/i?*+SME 

K (mag,ESO) 8.999 ± 0.017 2MASS C 

M K (mag,ESO) 2.36 ± 0.12 YY+a/fl*+SME 

Age (Gyr) 5.0+j^ YY+a/R*+SME 

Distance (pc) . . 214 ± 12 YY+a / ii* +SME 



a SME = "Spectroscopy Made Easy" package for 
analysis of high-resolution spectra Valenti & Piskunov 
(1996). These parameters depend primarily on SME, 
with a small dependence on the iterative analysis incor- 
porating the isochrone search and global modeling of 
the data, as described in the text. 

b SME+Claret = Based on SME analysis and tables 
from Claret (2004). 

c YY+a/iJ*+SME = Yi et al. (2001) isochrones, a/R* 
luminosity indicator, and SME results. 

d Based on the TASS catalogue (Droege et al. 2006). 

°Based on the transformations from Carpenter (2001). 

This section presents a simultaneous fitting of the 
HATNet photometry, the follow-up light curves, and the 
RV observations, which we refer to as "global" modeling. 
It incorporates not only a physical model of the system, 
but also a description of systematic (instrumental) vari- 
ations. 

Our model for the follow-up light curves used analytic 
formulae based on Mandel & Agol (2002) for the eclipse 
of a star by a planet, where the stellar flux is described 
by quadratic limb-darkening. The limb darkening coef- 
ficients were taken from the tables by Claret (2004) for 
the i band, corresponding to the stellar properties de- 
termined from the SME analysis (§ 3.1). The transit 
shape was parametrized by the normalized planetary ra- 
dius p = Rp/R*, the square of the impact parameter 
b 2 , and the reciprocal of the half duration of the tran- 
sit (/i?*. We chose these parameters because of their 
simple geometric meanings and the fact that they show 
negligible correlations (see e.g. Carter et al. 2008). 

Our model for the HATNet data was the simplified 
"P1P3" version of the Mandel & Agol (2002) analytic 
functions (an expansion by Legendre polynomials), for 
the reasons described in Bakos et al. (2009). The depth 
of the HATNet transits was adjusted independently in 
the fit (the depth was Bmst "blending factor" times that 
of the follow-up data) to allow for possible contamination 
by nearby stars in the under-sampled images of HATNet. 

As indicated earlier, initial modeling of the RV ob- 
servations showed deviations from a Kcplcrian fit highly 
suggestive of a second body in the system with a much 
longer period than the transiting planet. Thus, in our 



global modeling, the RV curve was parametrized by the 
combination of an eccentric Kcplcrian orbit for the inner 
planet with semi-amplitude K, and Lagrangian orbital 
elements (k, h) = ex (cosu, sino>), plus an eccentric Ke- 
plerian orbit for the outer object with K 2l fc 2 and h 2 , 
and a systemic RV zero-point 7. Throughout this paper 
the subscripts "1" and "2" will refer to HAT-P-13b and 
HAT-P-13c, respectively. If the subscript is omitted, we 
refer to HAT-P-13b. 

In the past, for single transiting planet scenarios we 
have assumed strict periodicity in the individual tran- 
sit times (Hartman ct al. 2009, and references therein), 
even when a drift in the RV measurements indicated an 
outer companion (HAT-P-llb, Bakos et al. 2009). Since 
the expected transit timing variations (TTV) for these 
planets were negligible compared to the error-bars of the 
transit center measurements, the strict periodicity was 
a reasonable hypothesis. Those data were characterized 
by two transit centers (Ta and Tb), and all interme- 
diate transits were interpolated using these two epochs 
and their corresponding N tr transit numbers. The model 
for the RV data component also implicitly contained its 
ephemeris information through Ta and Tb, and thus was 
coupled with the photometry data in the time domain. 

For HAT-P-13b, however, the disturbing force by the 
outer planet HAT-P-13c is expected to be not insignif- 
icant, because the RV semi-amplitude of the host star 
(<~ 0.5 kms -1 ) indicates that HAT-P-13c is massive, and 
the relatively short period and eccentric orbit (see later) 
indicate that it moves in relatively close to HAT-P-13b. 
Thus, the assumption of strict periodicity for HAT-P- 
13b is not precisely correct. While we performed many 
variations of the global modeling, in our finally adopted 
physical model we assume strict periodicity only for the 
HATNet data, where the timing error on individual tran- 
sits can be ~ 1000 seconds. In this final model we al- 
low for departure from such a periodicity for the indi- 
vidual transit times for the seven follow-up photometry 
observations. In practice, we assigned the transit num- 
ber N tr — to the first complete high quality follow- 
up light curve gathered on 2008 November 8/9 (MST). 
The HATNet data-set was characterized by T Cj _37o and 
T c , -309, covering all transit events observed by HATNet 
(here the c subscript denotes "center" for the transits of 
HAT-P-13b). The transit follow-up observations were 
characterized by their respective times of transit cen- 
ter: T^_68j T Cj _i, T Ct o, T c i, T c ,24, ?c,35: 7c, 62- Initial 
estimates of the T c ^ values yielded an initial epoch T c 
and period P by linear fitting weighted by the respec- 
tive error-bars of the transit centers. The model for the 
RV data component of the inner planet contained the 
ephemeris information through the T c and P variables, 
i.e. it was coupled with the transit data. The global mod- 
eling was done in an iterative way. After an initial fit to 
the transit centers (and other parameters; see later), T c 
and P were refined, and the fit was repeated. 

The time dependence of the RV of the outer planet 
was described via its hypothetical transit time T 2c (as 
if its orbit were edge on), and its period P%. The time 
of periastron passage T 2per i can be equivalently used in 
place of time of conjunction T 2c . 

Altogether, the 21 parameters describing the physi- 
cal model were r c ,_ 370 , T Ci _ 309 (HATNet), T Ci _ 68 , T c _i, 
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T c ,o, T Ci i, T c , 24 , T Ci35 , T c , 62 (FLWO 1.2m), b 2 , 
C/R*, K, 7, fc = ecosw, /i = esinw, i^, &2, /12, and 
T2c- Two additional auxiliary parameters were the in- 
strumental blend factor B inst of HATNct, and the HAT- 
Net out-of-transit magnitude, M 0j HATNet- 

-0.01 p 1 1 1 1 1 1 1 n 
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Time from transit center (days) 

Fig. 5. — Part of the global modeling described in § 3 are correc- 
tions for systematic variations of the light curves via simultaneous 
fitting with the physical model of the transit. The figure shows the 
resulting EPD- and TFA-corrected light curves for all 7 follow-up 
events (small gray points), and the merged and binned light curve 
(bin-size 0.005 days). 

We extended our physical model with an instrumen- 
tal model that describes the systematic (non-physical) 
variations (such as trends) of the follow-up data. This 
was done in a similar fashion to the analysis presented 
in Bakos et al. (2009). The HATNet photometry had 
already been EPD- and TFA-corrected before the global 
modeling, so we only modeled systematic effects in the 
follow-up light curves. We chose the "ELTG" method, 
i.e. EPD was performed in "local" mode with EPD co- 
efficients defined for each night, whereas TFA was per- 
formed in "global mode" , with the same coefficients de- 
scribing the optimal weights for the selected template 
stars in the field. The five EPD parameters include 
the hour angle (characterizing a monotonic trend that 
changes linearly over a night), the square of the hour 
angle, and three stellar profile parameters (equivalent to 
FWHM, elongation, and position angle). The exact func- 
tional form of the above parameters contained 6 coeffi- 
cients, including the auxiliary out-of-transit magnitude 
of the individual events. The EPD parameters were in- 
dependent for all 7 nights, implying 42 additional coef- 
ficients in the global fit. For the global TFA analysis 
we chose 19 template stars that had good quality mea- 
surements for all nights and on all frames, implying 19 
additional parameters in the fit. Thus, the total number 
of fitted parameters was 21 (physical model) + 2 (auxil- 
iary) + 42 (local EPD) + 19 (TFA) = 84. 

The joint fit was performed as described in Bakos et 
al. (2009), by minimizing x 2 in the parameter space us- 
ing a hybrid algorithm, combining the downhill simplex 
method (AMOEBA, see Press et al. 1992) with the clas- 
sical linear least-squares algorithm. We used the par- 
tial derivatives of the model functions as derived by Pal 
(2008). Uncertainties in the parameters were derived us- 
ing the Markov Chain Monte-Carlo method (MCMC, see 
Ford 2006). Since the eccentricity of the inner system 
appeared as marginally significant (k — —0.016 ± 0.008, 
h = 0.008 ± 0.012, implying e = 0.021 ± 0.009), and also 



because the physical model dictates that in the presence 
of a massive outer companion, the inner planet could 
maintain a non-zero eccentricity, we did not fix k and h to 
zero. The best fit results for the relevant physical param- 
eters are summarized in Tab. 4 and Tab. 5. The individ- 
ual transit centers for the FLWO 1.2 m data are given in 
§ 4.2. Other parameters fitted but not listed in the table 
were: T c _ 370 = 2453700.9182 ± 0.0068 (BJD), T c _ 309 = 
2453878.8216 ± 0.0093 (BJD), B mst = 0.82 ± 0.06, and 
M ,HATNot = 9.7966 ± 0.0001 (I band). The fit to the 
HATNet photometry data was shown earlier in Fig. 1, 
the orbital fit to the RV data is shown in Fig. 2, and 
the fit to the FLWO 1.2 m data is displayed in Fig. 3. 
The low rms of 3.4 ms -1 around the orbital fit is due in 
part to the use of an iodine free Keck/HIRES template 
that was acquired with higher spectral resolution and 
higher S/N than usual. Note that the low rms implies 
the absence of any additional interior planets other than 
HAT-P-13b, consistent with expectations given that the 
massive and eccentric outer planet HAT-P-13c dynam- 
ically forbids such planets (see e.g. Wittenmyer et al. 
2007). The stellar jitter required to reconcile the reduced 
X 2 with 1.0 for the RV data is 3.0ms _1 . The low jitter is 
also consistent with HAT-P-13 being a chromosphcrically 
quiet star (based on log R' HK and S) . 

The planetary parameters and their uncertainties can 
be derived by the direct combination of the a posteri- 
ori distributions of the light curve, radial velocity and 
stellar parameters. We find that the mass of the inner 
planet is M p = 0.853^;^ Mj, the radius is R p = 1.281± 
0.079 R 3 , and its density is p p = 0.498^;^ gem -3 . The 
final planetary parameters are summarized at the bottom 
of Table 4. The simultaneous EPD- and TFA-corrected 
light curve of all photometry follow-up events is shown 
in Fig. 5. 

The outer planet, HAT-P-13c, appears to be very mas- 
sive, with mi sin i 2 = 15.2±1.0 Mj, and orbits the star in 
a highly eccentric orbit with e 2 — 0.691 ±0.018. The ori- 
entation of the orbit (o>2 = 176.7 ±0.5°) is such that our 
line of sight is almost along the minor axis, coincidentally 
as in the HAT-P-2b system (Bakos et al. 2007). We note 
that the periastron passage of HAT-P-13c has not been 
well monitored, and thus the RV fit of the orbit suffers 
from a strong correlation between the quantities K 2 , e 2 
and 7, leading to correlated m 2 sinz and e 2 (Fig. 6). 



4. DISCUSSION 

We present the discovery of HAT-P-13, the first de- 
tected multi-planet system with a transiting planet. The 
inner transiting planet, HAT-P-13b, is an inflated "hot 
Jupiter" in a nearly circular orbit. The outer planet, 
HAT-P-13c, is both extremely massive and highly eccen- 
tric, and orbits in a P > 1 yr orbit . With an iron 
abundance of [Fe/H] = +0.41 ±0.08, the host star is also 
remarkable. As we describe below, this extraordinary 
system is a rich dynamical laboratory, the first to have 
an accurate clock (HAT-P-13b) and a known perturbing 
force (HAT-P-13c). The inner planet will help refine the 
orbital configuration (and thus true mass) of the outer 
planet through transit timing variations (TTVs). Con- 
versely, the outer planet, through its known perturba- 



8 



Bakos et al. 




TABLE 4 

Orbital and planetary parameters for HAT-P-13b 



m 2 sin(i) 

Fig. 6. — The orbital parameters, and as a consequence, the 
mass (rrt2 sin i) and eccentricity ei of the outer planet are strongly 
correlated (see § 3). Gray dots represent the results of the MCMC 
runs (20,000 trials). The final solution is marked with an open 
circle. The bimodal distribution in m2Sini is due to the similar 
distribution of the stellar mass (see § 3.1). 

tion, may constrain structural parameters and the tidal 
dissipation rate of the inner planet (Batygin, Bodcn- 
heimer, & Laughlin 2009), in addition to all the informa- 
tion that can be gleaned from the transits of HAT-P-13b. 

4.1. The planet HAT-P-13b 

The only other known planet with a host-star as 
metal rich as HAT-P-13 ([Fe/H]=+0.41 ± 0.08) is XO- 
2b ([Fe/H] = +0.45, Burke et al. 2007). XO-2b, how- 
ever, has much smaller mass (0.56 Mj) and smaller radius 
(0.98 Rj) than HAT-P-13b (0.85 Mj, 1.28 Rj). Other 
planets similar to HAT-P-13b include HAT-P-9b (M p = 
0.78 Mj, R. p = lARj, Shporer et al. 2008), and XO-lb 
(M p = 0.92 Mj, R p = 1.21 Rj, McCullough et al. 2006). 

When compared to theoretical models, HAT-P-13b is 
a clearly inflated planet. The Baraffe et al. (2008) mod- 
els with solar insolation (at 0.045 AU) are consistent 
with the observed mass and radius of HAT-P-13b ei- 
ther with overall metal content Z=0.02 and a very young 
age of 0.05-0.1 Gyr, or with metal content Z=0.1 and 
age 0.01-0.05 Gyr. 12 Given the fact, however, that the 
host star is metal rich, and fairly old (5.0to'.7 Gyr), it is 
unlikely that HAT-P-13b is newly formed (50Myr) and 
very metal poor. Comparison with Fortney et al. (2007) 
leads to similar conclusions. Using these models, HAT- 
P-13b is broadly consistent only with a 300 Myr planet 
at 0.02 AU solar distance and core-mass up to 25 M e , 
or a lGyr core-less (pure H/He) planet. If HAT-P-13b 
has a significant rocky core, consistent with expectation 
given the high metallicity of HAT-P-13, it must somehow 
be inflated beyond models calculated for such a planet 
with age and insolation suggested by the data. Nu- 
merous explanations have been brought up in the past 
to explain the inflated radii of certain extrasolar plan- 
ets (Miller, Fortney, & Jackson 2009; Fabrycky, John- 
son, & Goodman 2007, and references therein). One 
among these has been the tidal heating due to non-zero 
eccentricity (Bodenheimer, Lin, & Mardling 2001). No 
perturbing companion has been found for the inflated 
planets HD 209458b and HAT-P-lb (for an overview, 
see Mardling 2007). The HAT-P-13 system may be the 



Parameter 



Value 



Light curve parameters, HAT-P-13b 

P (days) 2.916260 ± 0.000010 

T c (BJD) 2454779.92979 ± 0.00038 

T14 (days) a 0.1345 ± 0.0017 

T12 = T 34 (days) a 0.0180 ± 0.0018 

a/fl* 5.84 ±0.26 

CAR* 17.07 ±0.16 

Rp/R* 0.0844 ± 0.0013 

b 2 4 47+ - 044 

b = a cos i/R+ 0.668l°-Q4g 

i (deg) 83.4 ± 0.6 

T peri (days) 2454780.64 ± 0.42 

RV parameters, as induced by HAT-P-13b 

K (ms^ 1 ) 106.1 ± 1.4 

k -0.016 ±0.008 

h 0.008 ±0.012 

e 0.021 ±0.009 

u> 181 ±45° 

Other RV parameters 

7 re! (ms- 1 ) -51.3 ±3.8 

Jitter (ms^ 1 ) 3.0 

Secondary eclipse parameters for HAT-P-13b 

T a (BJD) 2454781.359 ±0.014 

T s 14 0.1345 ±0.0069 

T 3 'i2 0.0180 ±0.0018 

Planetary parameters for HAT-P-13b 

M P (Mj) 0.853t°-°g 

Rp (Rj) 1.281 ±0.079 

C(M P ,R P ) b 0.56 

Pv (gcm-3) 0.4981°-^ 

a (au) ".,i.;27 

logc/ p (cgs) 3.11 ± 0.05 

T eq (K) 1653 ± 45 

6 c 0.046 ± 0.003 

Fp e r (ergs" 1 cm" 2 ) d (1.76 ± 0.20) ■ 10 13 

F ap (ergs" 1 cm- 2 ) c (1.62 ± 0.18) ■ 10 13 

(F) (ergs- 1 cm" 2 ) (1.69 ± 0.18) ■ 10 13 

a T14: total transit duration, time between first to last contact; 
T12 = T34: ingress/egress time, time between first and second, 
or third and fourth contact. 
b Correlation coefficient between the planetary mass M p and 



radius Rp 

c The Safronov number is given by © = 7j(V c£ 
(a/Rp)(M p /M*) (see Hansen & Barman 2007). 
d Incoming flux per unit surface area. 



:/V orb ) 2 



12 The equivalent semi-major axis, at which HAT-P-13b would 
receive the same amount of insolation when orbiting our Sun, is 

equiv 



a, 



: 0.0285 ± 0.0016 AU. 



first, where the inflated radius can be explained by the 
non-zero eccentricity of HAT-P-13b, excited by the outer 
HAT-P-13c planet orbiting on a highly eccentric orbit. 
We note that while the eccentricity of HAT-P-13b is non- 
zero only at the 2-a level (because k is non-zero at 2-<r), 
its pericenter is aligned (to within 4° ± 40°) with that 
of the outer planet HAT-P-13c. Additional RV measure- 
ments and/or space-based timing of a secondary eclipse 
are necessary for a more accurate determination of the or- 
bit. Nevertheless, if the apsidal lines of HAT-P-13b and 
HAT-P-13c are indeed aligned, and, if the configuration 
is co-planar, then the system is in a tidal fix-point config- 
uration, as recently noted by Batygin, Bodenheimer, & 
Laughlin (2009). This configuration imposes constraints 
on the structure of the inner planet HAT-P-13b. In par- 
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TABLE 5 

Orbital and planetary parameters for HAT-P-13C 



Parameter 



Value 



RV parameters, as induced by HAT-P-13c 

P 2 (days) 428.5 ± 3.0 

T 2c a (BJD) 2454870.4 ± 1.8 

K 2 (ms" 1 ) 502 ± 37 

k 2 -0.690 ±0.018 

h 2 0.039 ±0.005 

e 2 0.691 ± 0.018 

u) 2 176.7 ±0.5° 

T 2 , per i (days) 2454890.05 ± 0.48 

Fictitious light curve parameters, HAT-P-13c b 

T 2 i4 c (days) 0.621 ± 0.029 

T 2 ] 12 = T: i4 (days) 0.0355 ± 0.0012 

Fictitious secondary eclipse parameters for HAT-P-13c a 

T 2a (BJD) 2454912.7 ±1.8 

T 2a 14 (days) 0.574 ± 0.028 

T 2s \i2 (days) 0.0355 ± 0.0012 

Planetary parameters for HAT-P-13c 

m 2 sini 2 (Mj) 15.2 ±1.0 

oa (AU) l-188l°' J 8 3 

T 2 , oq (K) 340 ± 9 

F 2 , per (ergs" 1 cm" 2 ) d (2.26 ± 0.35) • 10 11 

F2,a P (ergs- 1 ™- 2 ) d (7.61 ±0.86) ■ 10 9 

(F 2 ) (ergs" 1 cm- 2 ) d (3.00 ± 0.33) ■ 10 10 

a T 2c would be the center of transit of HAT-P-13c, if its 
(unknown) inclination was 90° . 

b Transits of HAT-P-13c have not been observed. The 
values are for guidance only, and assume zero impact pa- 
rameter. 

c T14: total transit duration, time between first to last 
contact, assuming zero impact parameter. T\ 2 = T34: 
ingress/egress time, time between first and second, or third 
and fourth contact. Note that these values are fictitious, 
and transits of HAT-P-13c have not been observed. 

d Incoming flux per unit surface area in periastron, apas- 
tron, and averaged over the orbit. 

ticular, Batygin, Bodcnhcimer, & Laughlin (2009) give 
limits on the tidal Love number k2b, core- mass and tidal 
energy dissipation rate Qb of HAT-P-13b. 

4.2. Transit timing variations 

It has long been known that multiple planets in the 
same planetary system perturb each other, and this may 
lead to detectable variations in the transit time of the 
transiting planet(s). Transit timing variations have been 
analytically described by Holman & Murray (2005) and 
Agol et al. (2005), and have been suggested as one of the 
most effective ways of detecting small mass perturbers of 
transiting planets. This has motivated extensive follow- 
up of known TEPs e.g. by the Transit Light Curve (TLC) 
project (Holman et al. 2006; Winn et al. 2007), looking 
for companions of TEPs. However, for HAT-P-13b there 
is observational (spectroscopic) evidence for an outer 
component (HAT-P-13c), and thus transit timing vari- 
ations of HAT-P-13b must be present. Transit timing 
variations can be used to constrain the mass and orbital 
elements of the perturbing planet, in our case those of 
HAT-P-13c. 

The global modeling described in § 3.3 treats the tran- 
sit centers of the FLWO 1.2 m telescope follow-up as in- 
dependent variables, i.e., it automatically provides the 
basis for a TTV analysis. Throughout this work, by TTV 



we refer to the time difference between the: observed tran- 
sit center (T c ^) and the calculated value based on a fixed 
epoch and period as given in Tab. 4, i.e. in the O — C 
sense. The resulting TTVs are shown in Fig. 7. We 
believe the error-bars to be realistic as they are the re- 
sult of a full MCMC analysis, where all parameters are 
varied (including the EPD and TFA parameters). As 
further support for this, the transits around N tr = 
(T Ct -i,T Ct o,T Ct i) show a standard deviation that is com- 
parable to the error-bars (and we can safely assume zero 
TTV within a ±2.9163 day time range). It is also appar- 
ent from the plot that the smallest error-bars correspond 
to the full transit observations (N tr = and 24). TTVs 
of the order of 100 seconds are seen from the best fit pe- 
riod. Given the large error-bars on our transit centers (of 
the order of 100 seconds), in part due to possible remain- 
ing systematics in the partial transit light curve events, 
we consider these departures suggestive only. 

Nevertheless, it is tempting to compare our results with 
analytic formulae presented by Agol et al. (2005) and 
Borkovits et al. (2003). The HAT-P-13(b,c) system cor- 
responds to case "ii" of Agol et al. (2005), i.e. with an 
exterior planet on an eccentric orbit having a much larger 
semi-major axis than the inner planet. Formulae for the 
general (non co-planar) case are given by Borkovits et 
al. (2003, Eq. 46). The TTV effect will depend on the 
following known parameters for the HAT-P-13 system: 
M*, M p , P, i p , P2, e2, L02 1 m2sin«2 and T2 tPer i (these 
are listed in Tab. 4 and Tab. 5). The TTV will also 
depend on the following unknown parameters: the true 
inclination of HAT-P-13c, 22, and the relative angle of 
the orbital normals projected onto the plane of the sky, 
D = i} 1 — fl 2 (see Fig. 1 of Borkovits et al. 2003). In ad- 
dition to the gravitational perturbation of HAT-P-13c on 
HAT-P-13b, the barycenter of the inner subsystem (com- 
posed of the host star HAT-P-13, and the inner planet 
HAT-P-13b) orbits about the three-body barycenter due 
to the massive HAT-P-13c companion. This leads to 
light-travel time effects in the transit times of HAT-P- 
13b (TTV1) that are of the same order of magnitude as 
the TTV effect due to the perturbation (TTVg). 

We have evaluated the analytic formulae including 
both the TTVg and TTV1 effects for cases with i 2 = 
83.4 ± 0.6° (corresponding to co-planar inner and outer 
orbits, and to M 2 = 15.2 Mj), and i 2 = 8.1° (an ad 
hoc value yielding an almost face-on orbit with M 2 = 
105 Mj), and D = 0° or D = 90°. These four analytic 
models are illustrated in Fig. 7. The bottom panel of 
Fig. 7 shows the TTV1 and TTVg effects separately for 
the i 2 = 83.4±0.6° and D = 0° case. The i 2 = 83.4±0.6° 
(M 2 = 15.2 Mj) cases give TTV variations of the order of 
15 seconds. The functional form of the i2 = 83.4 ± 0.6° 
and D = 0° case appears to follow the observational 
values, albeit with much smaller amplitude. Curiously, 
for this case the TTV1 effect cancels the TTVg effect 
to some extent (bottom panel of Fig. 7). Increasing the 
mass of the outer companion by decreasing i 2 at constant 
7712 sini 2 does increase the TTV amplitude up to 100 sec- 
onds at 22 = 8.1°, but the functional form changes and 
no longer resembles the trend seen in the observational 
data. 

In conclusion, while it is premature to fit the cur- 
rent data-set with analytic models because of the small 
number of data-points (7) and the large error-bars (<~ 
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TABLE 6 
Transit timing variations of 
HAT-P-13. 



N tr 


T 
(BJD) 


cr-p 
(sec) 


O-C 

(sec) 


-68 


2454581.62406 


105.7 


-10.9 


-1 


2454777.01287 


86.7 


-58.5 





2454779.92953 


54.3 


-24.2 


1 


2454782.84357 


133.6 


-215.6 


24 


2454849.92062 


65.1 


50.7 


35 


2454882.00041 


129.5 


132.2 


62 


2454960.73968 


154.1 


156.1 



100 sec), these data are not inconsistent with the pres- 
ence of TTVs, and will prove useful in later analyses. Fu- 
ture measurements of full transit events with high accu- 
racy in principle can determine both i 2 and D, i.e., HAT- 
P-13c may become an RV-based detection with known 
orbital orientation and a true mass (even if it does not 
transit). 



i2=ii,n r n ? =o 

- : 2=M.<> .-LlU'JC 






i2=8.1.sVsii=° 

i2=8.1.£2 2 -Q 1 =90 































N. r transit number 



TTVg + TTVI 

« -TKF 




N lr transit number 

Fig. 7. — (Top) Transit times of the individual transit events 
of HAT-P-13b. The filled circles correspond to the global analysis 
described in § 4.2. The large filled square, the open circle and the 
open square correspond to the apastron, conjunction and pcrias- 
tron of the perturber HAT-P-13c. Overlaid are analytic models for 
the transit timing variation due to gravitational effects (TTVg) and 
light-travel time effects (TTVI) for four scenarios: i) the inclination 
of HAT-P-13c is i 2 = 83.1° = h, ii) i 2 = 8.1° (at fixed m 2 sini 2 , 
corresponding to a large mass for HAT-P-13c), iii) i 2 = 83.1° and 
the mutual inclination of the orbital normals in the plane of the 
sky is D = Hi - U 2 = 90° (see Fig. 1 of Borkovits et al. (2003)), or 
iv) i 2 = 8.1° and D = 90°. (Bottom) The selected case "i" from 
above with zoomed-in vertical scale, showing the TTVg (dashed 
line) and TTVI (dotted line) effects separately, and their net effect 
(solid line). 



4.3. The planet HAT-P-13c 

The probability of HAT-P-13c transiting the host star, 
as seen from the Earth, depends on R+, R 2p , e 2 , uj 2 , 
and a 2 (Kane & von Braun 2009). We evaluated the 
transit probability in a Monte Carlo fashion, as part of 
the global modeling, resulting in P 2 . tr = 0.0130 ± 0.0008 



if R 2 , p = lRj (P 2 . tr = 0.0122 ± 0.0007 if R 2 . p -> 0). 
This derivation assumes an isotropic inclination distri- 
bution. However, dynamical constraints, such as precise 
measurements of TTV effects, or analysis of orbital sta- 
bility, may further limit the allowed inclination range, 
and thus increase or decrease the chance of transits. 

Unfortunately, our HATNet and FLWO 1.2 m datasets 
do not cover any time interval around the expected tran- 
sit times of HAT-P-13c, thus we can not prove or refute 
the existence of such transits. Nevertheless, it is an inter- 
esting thought experiment to characterize the putative 
transits of HAT-P-13c, should they occur. HAT-P-13c 
would be the longest period transiting planet discovered 
to date, with a semi-major axis of over 1 AU, and a pe- 
riod of ~ 428 days, about 4 times longer than the current 
record holder HD 80606 (Naef et al. 2001; Fossey, Wald- 
mann, & Kipping 2009; Moutou et al. 2009). With a 
true mass of 15.2 Mj, we have good reasons to believe 
that its radius would be around based on heavy- 

mass transiting planets like HAT-P-2b (8.84 Mj, Bakos 
et al. 2007), XO-3b (11.79 Mj, Johns-Krull et al. 2008), 
Corot-3b (21.66 Mj, Deleuil et al. 2008), all having radii 
around 1 Rj . If the transits are full (i.e. not grazing) , 
then the transit depth would be similar to that of the 
inner planet HAT-P-13b. The duration of the transit 
could be up to 14.9 hours. Follow-up observations of such 
transits would require cither a world-wide effort, or un- 
interrupted space-based observations. The next transit 
center, according to the present analysis, will occur at 
2010 April 12 9am UTC. Since we are looking at the or- 
bit of HAT-P-13c along the semi-latus rectum (parallel 
to the minor- axis), the chance for an occultation of the 
planet by the star has a very similar chance of occurancc 
as the primary eclipse. Observations of the secondary 
eclipse could greatly decrease the error on e 2 and uj 2 . 

As regards to the nature of HAT-P-13c, even at its 
minimal mass (15.2 ± 1.0 Mj) it is the 10th most massive 
planet out of 327 planets listed on the Extrasolar Planet 
Encyclopedia as of 2009 July. Curiously, the recently 
announced Doppler-detection HD 126614b (Howard et 
al. 2009) appears to have similar orbital characteristics 
to HAT-P-13c. Both are Jovian planets in P > lyr, 
e = 0.7 orbits around metal-rich ([Fe/H] « 0.5) stars. 
As described earlier in § 4.2, we have good hopes that 
in the near future precise TTV measurements of the in- 
ner planet HAT-P-13b, will constrain the orbital incli- 
nation and thus the true mass of HAT-P-13c. Further, 
such TTV variations can also constrain the mutual in- 
clination of HAT-P-13b and HAT-P-13c. Measuring the 
sky-projected angle between the stellar spin axis and the 
orbital normal of HAT-P-13b (the inner planet) via the 
Rossitcr-McLaughlin effect will shed light on the migra- 
tion history of HAT-P-13b, and by inference the scatter- 
ing history between HAT-P-13b and HAT-P-13c. 
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